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Abstract 

We present a numerical method to compute the optimal maintenance time for a complex dynamic 
system applied to an example of maintenance of a metallic structure subject to corrosion. An arbitrarily 
early intervention may be uselessly costly, but a late one may lead to a partial/complete failure of 
the system, which has to be avoided. One must therefore find a balance between these too simple 
maintenance policies. To achieve this aim, we model the system by a stochastic hybrid process. The 
maintenance problem thus corresponds to an optimal stopping problem. We propose a numerical method 
to solve the optimal stopping problem and optimize the maintenance time for this kind of processes. 

Index Terms 

Dynamic reliability, predictive maintenance, Piece-wise-deterministic Markov processes, optimal 
stopping times, optimization of maintenance. 

I. Introduction 

A complex system is inherently sensitive to failures of its components. We must therefore 
determine maintenance policies in order to maintain an acceptable operating condition. The 
optimization of maintenance is a very important problem in the analysis of complex systems. It 
determines when maintenance tasks should be performed on the system. These intervention dates 
should be chosen to optimize a cost function, that is to say, maximize a performance function 
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or, similarly, to minimize a loss function. Moreover, this optimization must take into account 
the random nature of failures and random evolution and dynamics of the system. Theoretical 
study of the optimization of maintenance is also a crucial step in the process of optimization of 
conception and study of the life service of the system before the first maintenance. 

We consider here an example of maintenance related to an aluminum metallic structure subject 
to corrosion. This example was provided by Astrium. It concerns a small structure within a 
strategic ballistic missile. The missile is stored successively in a workshop, in a nuclear submarine 
missile launcher in operation or in the submarine in dry-dock. These various environments are 
more or less corrosive and the structure is inspected with a given periodicity. It is made to have 
potentially large storage durations. The requirement for security is very strong. The mechanical 
stress exerted on the structure depends in part on its thickness. A loss of thickness will cause an 
over-constraint and therefore increase a risk of rupture. It is thus crucial to control the evolution 
of the thickness of the structure over time, and to intervene before the failure. 

The only maintenance operation we consider here is the complete replacement of the structure. 
We do not allow partial repairs. Mathematically, this problem of preventive maintenance corre- 
sponds to a stochastic optimal stopping problem as explained by example in the book of Aven 
and Jensen [1 J. It is a difficult problem, because on the one hand, the structure spends random 
times in each environment, and on the other hand, the corrosiveness of each environment is also 
supposed to be random within a given range. In addition, we search for an optimal maintenance 
date adapted to the particular history of each structure, and not an average one. We also want to 
be able to update the predicted maintenance date given the past history of the corrosion process. 

To solve this maintenance problem, we propose to model this system by a piecewise-deter- 
ministic Markov process (PDMP). PDMP's are a class of stochastic hybrid processes that have 
been introduced by Davis (31 in the 80's. These processes have two components: a Euclidean 
component that represents the physical system (e.g. temperature, pressure, thickness loss) and a 
discrete component that describes its regime of operation and/or its environment. Starting from 
a state x and mode m at the initial time, the process follows a deterministic trajectory given 
by the laws of physics until a jump time that can be either random (e.g. it corresponds to a 
component failure or a change of environment) or deterministic (when a magnitude reaches a 



February 7, 2011 



DRAFT 



3 



certain physical threshold, for example the pressure reaches a critical value that triggers a valve). 
The process restarts from a new state and a new mode of operation, and so on. This defines 
a Markov process. Such processes can naturally take into account the dynamic and uncertain 
aspects of the evolution of the system. A subclass of these processes has been introduced by 
Devooght [5] for an application in the nuclear field. The general model has been introduced in 
dynamic reliability by Dutuit and Dufour ||6|. 

The theoretical problem of optimal stopping for PDMP's is well understood, see e.g. Gugerli 
0. However, there are surprisingly few works in the literature presenting practical algorithms to 
compute the optimal cost and optimal stopping time. To our best knowledge only Costa and Davis 
13 have presented an algorithm for calculating these quantities for PDMP's. Yet, as illustrated 
above, it is crucial to have an efficient numerical tool to compute the optimal maintenance time 
in practical cases. The purpose of this paper is to adapt the general algorithm recently proposed 
by the authors in JH to this special case of maintenance and show its high practical power. 
More precisely, we present a method to compute the optimal cost as well as a quasi optimal 
stopping rule, that is the date when the maintenance should be performed. As a byproduct of our 
procedure, we also obtain the distribution of the optimal maintenance dates and can compute 
dates such that the probability to perform a maintenance before this date is below a prescribed 
threshold. 

The remainder of this paper is organized as follows. In section |TTJ. we present the example 
of corrosion of the metallic structure that we are interested in with more details as well as the 



framework of PDMP's. In section III we briefly recall the formulation of the optimal stopping 



problem for PDMP's and its theoretical solution. In section IV, we detail the four main steps of 



algorithm. In section [V] we present the numerical results obtained on the example of corrosion. 



Finally, in section VI, we present a conclusion and perspectives. 



II. Modeling 

Throughout this paper, our approach will be illustrated on an example of maintenance of a 
metallic structure subject to corrosion. This example was proposed by Astrium. As explained 
in the introduction, it is a small homogeneous aluminum structure within a strategic ballistic 
missile. The missile is stored for potentially long times in more or less corrosive environments. 
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The mechanical stress exerted on the structure depends in part on its thickness. A loss of thickness 
will cause an over-constraint and therefore increase a risk of rupture. It is thus crucial to control 
the evolution of the thickness of the structure over time, and to intervene before the failure. 

Let us describe more precisely the usage profile of the missile. Its is stored successively in 
three different environments, the workshop, the submarine in operation and the submarine in 
dry-dock. This is because the structure must be equipped and used in a given order. Then it goes 
back to the workshop and so on. The missile stays in each environment during a random duration 
with exponential distribution. Its parameter depends on the environment. At the beginning of 
its service time, the structure is treated against corrosion. The period of effectiveness of this 
protection is also random, with a Weibull distribution. The thickness loss only begins when 
this initial protection is gone. The degradation law for the thickness loss then depends on the 
environment through two parameters, a deterministic transition period and a random corrosion 
rate uniformly distributed within a given range. Typically, the workshop and dry-dock are the 
more corrosive environments. The randomness of the corrosion rate accounts for small variations 
and uncertainties in the corrosiveness of each environment. 

We model this degradation process by a 3-dimensional PDMP (X t ) with 3 modes correspond- 
ing to the three different environment. Before giving the detailed parameters of this process, we 
shortly present general PDMP's. 

A. Definition of piecewise-deterministic Markov processes 

Piecewise-deterministic Markov processes (PDMP's) are a general class of hybrid processes. 
Let M be the finite set of the possible modes of the system. In our example, the modes correspond 
to the various environments. For all mode m in M, let E m an open subset in R d . A PDMP is 
defined from three local characteristics ($, A, Q) where 

• the flow $ : M x R d x R ->■ R d is continuous and for all s,t>0, one has $(•, •, t + s) = 
$($(-,-, s),t). It describes the deterministic trajectory of the process between jumps. For 
all (m, x) in M x E m , we set 

t*(m,x) = inf{t > : $(m,x,t) G dE m }, 

the time to reach the boundary of the domain starting from x in mode m. 
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• the jump intensity A characterizes the frequency of jumps. For all (m, x) in M x E m , and 

t < t*(m, x), we set 



• the Markov kernel Q represents the transition measure of the process and allows to select 

the new location after each jump. 
The trajectory X t = (m t ,x t ) of the process can then be defined iteratively. We start with an 
initial point X = (k ,y ) with k e M and y e E ko . The first jump time Ti is determined by 



On the interval [0,Ti), the process follows the deterministic trajectory m t = k and x t = 
$(/c , y , t). At the random time Ti, a jump occurs. Note that a jump can be either a discontinuity 
in the Euclidean variable x t or a change of mode. The process restarts at a new mode and/or 
position X Tl = (kx,yi), according to distribution Q ko (Q(k , y , Ti), •). We then select in a similar 
way an inter jump time T 2 — Ti, and in the interval [Ti, T 2 ) the process follows the path m t = k\ 
and xt = &(ki,yi,t — T\). Thereby, iteratively, a PDMP is constructed, see Figure [T] for an 
illustration. Let Z = X , and for n > 1, Z n = X Tn , location and mode of the process after 



Fig. 1. An exemple of path for a PDMP until the second jump. The first jump is random. The second jump is deterministic 
because the process has reached the boundary of the domain. 

each jump. Let So = 0, S\ = Ti and for n > 2, S n — T n —T n _i the inter-jump times between two 
consecutive jumps, then (Z n , S n ) is a Markov chain, which is the only source of randomness of 
the PDMP and contains all information on its random part. Indeed, if one knows the jump times 
and the positions after each jump, we can reconstruct the deterministic part of the trajectory 
between jumps. It is a very important property of PDMP's that is at the basis of our numerical 
procedure. 
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B. Example of corrosion of metallic structure 

We can now turn back to our example of corrosion of structure and give the characteristics 
of the PDMP modeling the thickness loss. The finite set of modes is M = {1,2,3}, where 
mode 1 corresponds to the workshop environment, mode 2 to the submarine in operation and 
mode 3 to the dry-dock. Although the thickness loss is a one-dimensional process, one needs 
a three dimensional PDMP to model its evolution, because it must also take into account all 
the sources of randomness, that is the duration of the initial protection and the corrosion rate in 
each environment. The corrosion process (X t ) is defined by: 



where m t is the environment at time t, d t is the thickness loss at time t, j t is the remainder of 
the initial protection at time t and p t is the corrosion rate of the current environment at time t. 

Originally, at time 0, one has X = (l,0,7 ,po)> which means that the missile is in the 
workshop and the structure has not started corroding yet. The original protection 7 is drawn 
according to a Weibull distribution function 



with a = 2.5 and (5 = 11800 hours -1 . The corrosion rate in the workshop is drawn according 
to a uniform distribution on [10~ 6 , 10~ 5 ] mm/hour. The time Ti spent in the workshop is drawn 
according to an exponential distribution with parameter Ai = 17520 hour -1 . At time t between 
time and time 7\, the remainder of the protection is simply j t = max{0, 70 — t}, p t is constant 
equal to p and the thickness loss d t is given by 



where % = 30000 hours. 

At time Ti, a jump occurs, which means there is a change of environment and a new corrosion 
rate is drawn for the new environment. The other two components of the process (X t ) modeling 
the remainder of the protection 7 t and the thickness loss d t naturally evolve continuously. 
Therefore, one has m Tl = 2, 77^ = if 70 < Ti, 7^ = 70 — Ti otherwise ; that is to say 
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that once the initial protection is gone, it has no effect any longer, p Tl is drawn according to 
a uniform distribution on [1CT 7 , 1CT 6 ] mm/hour. The process continues to evolve in the same 
way until the next change of environment occurring at time T 2 . Between T\ and T 2 , just replace 
Po by Pn, 7o by 7^, Vi by 772 = 200000 hours and t by t — T\ in equation ([TJ). The process 
visits successively the 3 environments always in the same order 1, 2 and 3 and then returns to 
the environment 1. . The time spent in the environment i is a random variable exponentially 
distributed with parameters Aj with Ai = 17520 hours -1 , A 2 = 131400 hours -1 and A 3 = 8760 
hours -1 . The thickness loss evolves continuously according to equation ([!]) with suitably changed 
parameters. The period of transition in the mode i is r\i with 771 = 30000 hours, 772 = 200000 hours 
and 773 = 40000 hours. The corrosion rate pi expressed in mm per hour is drawn at each change 
of environments. In environments 1 and 3, it follows a uniform distribution on [10~ 6 , 10~ 5 ] and 
in environment 2, it follows a uniform distribution on [10 -7 ,10 -6 ]. Figure [2] shows examples 
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(a) One trajectory (b) 100 trajectories 

Fig. 2. Examples of trajectories of thickness loss over time. 

of simulated trajectories of the thickness loss. The slope changes correspond to changes of 
environment. The observed dispersion is characteristic of the random nature of the phenomenon. 
Note that the various physical parameters were given by Astrium and will not be discussed here. 

The missile is inspected and the thickness loss of the structure under study is measured at each 
change of environment. Note that the structure is small enough for only one measurement point 
to be significant. The structure is considered unusable if the loss of thickness reaches 0.2mm. 
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The optimal maintenance time must therefore occur before reaching this critical threshold, 
which could cause the collapse of the structure, but not too soon which would be unnecessarily 
expensive. It should also only use the available measurements of the thickness loss. 



We now briefly formulate the general mathematical problem of optimal stopping corresponding 
to our maintenance problem. Let z = (k , y ) be the starting point of the PDMP (X t ). Let A4n 
be the set of all stopping times T for the natural filtration of the PDMP (X t ) satisfying T < T N 
that is to say that the intervention takes place before the iVth jump of process. The iVth jump 
represents the horizon of our maintenance problem, that is to say that we impose to intervene 
no later than iVth change of environment. The choice of iV is discussed below. Let g be the cost 
function to optimize. Here, g is a reward function that we want to maximize. The optimization 
problem to solve is the following 



The function v is called the value function of the problem and represents the maximum perfor- 
mance that can be achieved. Solving the optimal stopping problem is firstly to calculate the value 
function, and secondly to find a stopping time r that achieves this maximum. This stopping time 
is important from the application point of view since it corresponds to the optimum time for 
maintenance. In general, such an optimal stopping time does not exist. We then define e-optimal 
stopping times as achieving optimal value minus e, i.e. v(z) — e. 

Under fairly weak regularity conditions, Gugerli has shown in [7J that the value function v 
can be calculated iteratively as follows. Let vn = g be the reward function, and we iterate an 
operator L backwards. The function v thus obtained is equal to the value function v. 



1v N = g, 
v k = L(v k+ i,g), 0<k<N-l. 
The operator L is a complex operator which involves a continuous maximization, conditional 
expectations and indicator functions, even if the cost function g is very regular. 

L{w,g)(z)= sup {E[w(Z 1 )l Sl<uM * {z) + g(^(z,u))l Sl > U At*(z)\Z = z]} 

u<t*(z) (2) 

WE [w{Z x )\Z Q = z]. 



III. Optimal stopping problem 



v(z) = sup E z [g{X T )\ . 



T&M N 
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However, we can see that this operator depends only on the discrete time Markov chain (Z n , S n ). 
Gugerli also proposes an iterative construction of e-optimal stopping times, which is a bit too 
tedious and technical to be described here, see for details. 

For our example of metallic structure, we choose an arbitrary reward function that depends 
only on the loss of thickness, since this is the critical factor to monitor. Note that we could take 
into account the other components of our process without any additional difficulty. The reward 
function is built to reflect the fact that beyond a loss of thickness of 0.2mm, the structure is 
unusable, so it is too late to perform maintenance. Conversely, if the thickness loss is small, 
such a maintenance is unnecessarily costly. We use a piecewise affine function g which values 
are given at the points in the table in Figure |3j As for the choice of the computational horizon 




Fig. 3. Graphical representation and definition of the cost function as a function of the thickness loss 

N, numerical simulations show that over 25 changes of environment, all trajectories exceed the 
critical threshold of 0.2mm. We will therefore set the time horizon to be the 25th jump (N = 25). 

IV. Numerical procedure 

It is natural to propose an iterative algorithm to calculate an approximation of the value 
function based on a discretization of the operator L defined in equation Q. This poses several 
problems, related to maximizing continuous functions, the presence of the indicator and the 
presence of conditional expectations. We nevertheless managed to overcome these three problems, 
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using the specific properties of PDMP's, and in particular the fact that the operator L depends 
only on the Markov chain (Z n , S n ). Our algorithm for calculating the value function is divided 
into three stages described below: a quantization of the Markov chain (Z n , S n ), a path-adapted 
time discretization between jumps, and finally a recursive computation of the value function v. 
Then, the calculation of quasi-optimal stopping time only uses comparisons of quantities already 
calculated in the approximation of the value function, which makes this technique particularly 
attractive, see flU for more mathematical details. 

A. Quantization 

The goal of the quantization step is to replace the continuous state space Markov chain (Z n , S n ) 
by a discrete state space chain (Z n ,S n ). The quantization algorithm is described in details in 
e-g- L8J flU ifTOl or |[TT1|. The principle is to obtain a finite grid adapted to the distribution of the 
random variable, rather than building an arbitrary regular grid. We discretize random variables 
rather than the state space, the idea is to put more points in the areas oh high density of the 
random variable. The quantization algorithm is based on Monte Carlo simulations combined 
with a stochastic gradient method. It provides N + 1 grids T n , < n < N of dimension d + 2, 
one for each couple (Z n , S n ), with K points in each grid. The algorithm also provide weights 
for the grid points and probability transition between two points of two consecutive grids. 

We note p n the projection to the nearest neighbor (for the Euclidean norm) from IR d+2 onto 
T n . The approximation of the Markov chain (Z n , S n ) is constructed as follows: 

(Z n , S n ) = Pni^Zn, S n ). 

Note that Z n and S n depend on both Z n and S n . The quantization theory ensures that the L 2 
norm of the distance between (Z n , S n ) and (Z n , S n ) tends to as the number of points K in 
the quantization grids tends to infinity, see IfTOl . 

It should be noted that when the dimension of Z is large, iV is large and we want to obtain 
grids with a large number K of points, the quantization algorithm can be time-consuming. 
However, we can make this grids calculation in advance and store them. They depend only 
on the distribution of the process, and not on the cost function. Figure [4] gives an example of 
quantization grid for the standard normal distribution in two dimensions. It illustrates that the 
quantization algorithm puts more points in areas of high density. 
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B. Time discretization 

We now wish to replace the continuous maximization of the operator L by a finite maxi- 
mization, that is to say that we must discretize the time intervals [0, £*(;?)] for each z in the 
quantization grids. For this, we choose a time step A < t*(z) (which may depend on z) and we 
construct the grids G(z) = {t\, ■ ■ ■ ,£ n (*)} defined by 

• n(z) is the integer part minus 1 of t*(z)/A, 

• for 1 < i < n(z), U = iA. 

We obtain grids that not only do not contain t*(z), but in addition, their maximum is strictly 
less than t*(z) — A, which is a crucial property to derive error bounds for our algorithm, see 
flU. Note also that we only need a finite number of grids G(z), corresponding to the z in the 
quantization grids (T n ) < n < N . Calculation of this time grids can still be made in advance. 
Another solution is to store only A and n(z) which are sufficient to reconstruct the grids. 

In practice, we choose a A that does not depend on z. To ensure that we have no empty grid, 
we first calculate the minimum of t*(z) on all grids of quantization, then we choose a A adapted 
to this value. 
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C. Approximate calculation of the value function 

We now have all the tools to provide an approximation of the operator L. For each 1 < n < N, 
and for all z in the quantization grid at time n — 1, we set 

L n (w,g)(z) = maxjs w(i„_i)l 5n<1<At . w +^($(i ft _i,u))l Sii > ttAt . w |i ft _ 1 = z } 

w(Z n )\Z n _i = 



u<G(z) 

VE 



Note that because we have different quantized approximations at each time step, we also have 
different discretizations of operator L at each time step. We then construct an approximation of 
the value function by backward iterations of the L n : 

vn = g, 

-i) — L n {v n , g){Z n . x ), l<n<N. 

Then we take v (Z ) = Vq(z) as an approximation of the value function v at the starting point z 
of the PDMP. It should be noted that the conditional expectations taken with respect to a process 
with discrete state space are actually finite weighted sums. 

Theorem 4.1: Under assumptions of Lipschitz regularity of the cost function g and local 
characteristics (<&, A, Q) of the PDMP, the approximation error in the calculation of the value 
function is 

\\v (z)-v (z)\\ 2 <Cy/EQ 

where C is an explicit constant which depends on the cost function and local characteristics of 
the PDMP, and EQ is the quantization error. 

Since the quantization error tends to when the number of points in the quantization grid 
increases, this result shows the convergence of our procedure. Here, the order of magnitude as 
the square root of the quantization error is due to the presence of indicator functions, which 
slow convergence because of their irregularity. To get around the fact that these functions are 
not continuous, we use the fact that the sets where they are actually discontinuous are of very 
low probability. The precise statement of this theorem and its proof can be found in J4J. 

D. Calculation of a quasi-optimal stopping time 

We have also implemented a method to compute an e-optimal stopping time. The discretization 
is much more complicated and subtle than that of operator L, because we need both to use the 
true Markov chain (Z n , S n ) and its quantized version (Z n , S„). The principle is as follows: 
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• At time 0, with the values Z = z and So = 0, we calculate a first date R\ which depends 
on Zq, So and on the value that has realized the maximum in the calculation of Li(v\,g). 

• We then allow the process to run normally until the time Ri A 7\, that is the minimum 
between this computed time R± and the first change of environment. If Ri comes first, it 
is the date of near-optimal maintenance, if 7\ comes first, we reset the calculation. 

• At time T\, with the values of Z\ and Si, we calculate the second date R 2 which depends 
on Z\ and Si and on the the value that has realized the maximum in the calculation of 

L2(h,g)- 

• We then allow the process to run normally until the time (Ti +_R 2 ) AT 2 , that is the minimum 
between the computed remaining time R 2 and the next change of environment. If 7\ + R 2 
comes first, it is the date of near-optimal maintenance, if T 2 comes first, we reset the 
calculation, and so on until the iVth jump time where maintenance will be performed if it 
has not occurred before. 

We have also proved the quality of this approximation by comparing the expectation of the cost 
function of the process stopped by the above strategy to the true value function. This result, its 
proof and the precise construction of our stopping time procedure can be found in |@]. 

This stopping strategy is interesting for several reasons. First, this is a real stopping time for 
the original PDMP which is a very strong result. Second, it requires no additional computation 
compared to those made to approximate the value function. This procedure can be easily 
performed in real time, and only requires an observation of the process at the times of change of 
environment, which is exactly the available inspection data for our metallic structure. Moreover, 
even if the original problem is an optimization on average, this stopping rule is path-wise and is 
updated when new data arrive on the history of the process at each change of environment. 
Finally, as our stopping procedure is of the form intervene at such date if no change of 
environment has occurred in the meantime, it allows in some measure to have maintenance 
scheduled in advance, In particular, our procedure ensures that there will be no need to perform 
maintenance before a given date, which is crucial for our example as a submarine in operation 
should not be stopped at short notice. 
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V. Numerical results 

We have implemented this procedure for the optimization of the maintenance of the metallic 
structure described in section |TTJ. With our choice of reward function, it is easy to see that the 
true value function at z = is 4, which is the maximum of the reward function g, and an 
optimal stopping time is the first moment when the loss reaches 0.18 mm thick (value where 
g reaches its maximum). This is because the cost function only depends on the thickness loss, 
which evolves continuously increasingly over time. However, our numerical procedure is valid 
for any sufficiently regular reward function, and we shall not use the knowledge of the true 
value function or optimal stopping time in our numerical procedure. Besides, we recall that the 
thickness loss is not measured continuously. 

While running the algorithm described in the previous section, e encountered an unexpected 
difficulty for the construction of the quantization grids. Indeed, the scales of the different variables 
of the problem are radically different: from about 1CT 6 for p to 10 5 for the average time spent 
in environment 2. This poses a problem in the classical quantization algorithm as searching 
the nearest neighbor and gradient calculations are done in Euclidean norm, regardless of the 
magnitudes of the components. Figure [5] illustrates this problem by presenting two examples of 




(a) Classical algorithm (b) Algorithm with weighted Euclidean norm 

Fig. 5. Quantization grids for a uniform distribution on [0, 1] x [0,5000] 



quantization grids for a uniform distribution on [0, 1] x [0, 5000]. The left image shows the result 
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obtained by the conventional algorithm, the right one is obtained by weighting the Euclidean norm 
to renormalize each variable on the same scale. It is clear from this example that the conventional 
method is not satisfactory, because the grid obtained is far from uniform. This defect is corrected 
by a renormalization of the variables. We therefore used a weighted Euclidean norm to quantify 
the Markov chain associated with our degradation process. 

Figure [6] shows some projections of the quantization grids with 2000 points that we obtained. 
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(a) Environment 2, time Ti 



(b) Environment 3, time T2 



(c) Environment 2, time T10 
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(d) Environment 1, time T15 



(e) Environment 2, time T±g 



(f) Environment 2,time T25 



Fig. 6. Quantization grids with 2000 points for the inter-jump time (abscissa) and the thickness loss (ordinate). The scale 
changes for each graph. 



The times are chosen in order to to illustrate the random and irregular nature of the grids, they 
are custom built to best approach the distribution of the degradation process. 

Figure [7] shows two examples of computation of the quasi optimal maintenance time on two 
specific simulated trajectories. The thick vertical line represents the moment provided by the 
algorithm to perform maintenance. The other vertical lines materialize the moments of change 
of environment, the horizontal dotted line the theoretical optimum. In both examples, we stop 
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Fig. 7. Examples of stopped trajectories with the optimal maintenance time calculated by the algorithm. 

at a value very close to the optimum value. In addition, the intervention did take place before 
the critical threshold of 0.2mm. 

We calculated an approximate value function v in two ways. The first one is the direct 
method obtained by the algorithm described above. The second one is obtained by Monte Carlo 
simulation using the quasi-optimal stopping time provided by our procedure. The numerical 
results we obtained are summarized in Table [1} We see as expected, that the greater the number of 
points in the quantization grid, the better our approximation becomes. Furthermore, the specific 
form of this cost function g indicates that at the threshold of 1, the intervention takes place 
between 0.15 and 0.2mm, and when the threshold increases, this range is narrowed. We can 
therefore state that our approximation is good even for low numbers of grid points. The last 
column of the table also shows the validity of our stopping rule. It should be noted here that 
this rule does not use the optimal stopping time stop at the first moment when the thickness 
loss reaches 0.18mm. The method we use is general and implementable even when the optimal 
stopping time is unknown or does not exist. 

Moreover, we can also construct a histogram (Figure [8]) of the values of our stopping time, that 
is to say, a histogram of the values of effective moments of maintenance. We can also estimate 
the probability that this moment is below certain thresholds. These results are interesting for 
Astrium in the design phase of the structure to optimize margins from the specifications and to 
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Number of points 
in the quantization 
grids 


Approximation of the 
value function by the 
direct algorithm 


Approximation of the value 
function by Monte Carlo with the 
cjua si -optimal stopping time 


1 n 


1 AS 


o qa 

U.yi 






1 SA 
1 .01 


i on 

1UU 


9 QA 


9 i n 


zuu 


1 HQ 




500 


3.39 


3.15 


1000 


3.56 


3.43 


2000 


3.70 


3.60 


5000 


3.82 


3.73 


8000 


3.86 


3.75 



TABLE I 

Numerical results for the calculation of the value function. 




(a) Histogram of 100000 values of the optimal maintenance (b) Quantiles. 

time expressed in years. 

Fig. 8. Distribution and quantiles of the quasi-optimal stopping time. 



consolidate the design margins available. Thus, we can justify that with a given probability no 
maintenance will be required before the termination date of the contract. 
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VI. Conclusion 

We have applied the numerical method described in [4] on a practical industrial example to 
approximate the value function of the optimal stopping problem and a quasi-optimal stopping 
time for a piecewise-deterministic Markov process, that is the quasi optimal maintenance date 
for our structure. The quantization method we propose can sometimes be costly in computing 
time, but has a very interesting property: it can be calculated off-line. Moreover it depends 
only on the evolutionary characteristics of the model, and not on the cost function chosen, or 
the actual trajectory of the specific process we want to monitor. The calculation of the optimal 
maintenance time is done in real time. This method is especially attractive as its application 
requires knowledge of the system state only at moments of change of environment and not 
in continous time. The optimal maintenance time is updated at the moments when the system 
switches to another environment and has the form intervene at such date if no change of mode 
takes place in the meantime, which allows to schedule maintenance services in advance. 

We have implemented this method on an example of optimization of the maintenance of a 
metallic structure subject to corrosion, and we obtained very satisfactory results, very close to 
theoretical values, despite the relatively large size of the problem. These results are interesting 
for Astrium in the design phase of the structure to maximize margins from the specifications 
and to consolidate the avaible dimensional margins. Thus, we propose tools to justify that with 
a given probability that no maintenance will be required before the end of the contract. 

The application that we have presented here is an example of maintenance as good as new 

of the system. The next step will be to allow only partial repair of the system. The problem 

will then be to find simultaneously the optimal times of maintenance and optimal repair levels. 

Mathematically, it is an impulse control problem, which complexity exceeds widely that of the 

optimal stopping. Here again, the problem is solved theoretically for PDMP, but there is no 

practical numerical method for these processes in the literature. We now work in this direction 

and we hope to be able to extend the results presented above. 
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